GW quasi-particle spectra from occupied states only 
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We introduce a method that allows for the calculation of quasi-particle spectra in the GW ap- 
proximation, yet avoiding any explicit reference to empty one-electron states. This is achieved by 
expressing the irreducible polarizability operator and the self-energy operator through a set of lin- 
ear response equations, which are solved using a Lanczos-chain algorithm. We first validate our 
approach by calculating the vertical ionization energies of the benzene molecule and then show its 
potential by addressing the spectrum of a large molecule such as free-base tetraphenylporphyrin. 
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In spite of the formidable success met over the past 
forty years in the simulation of materials, based on 
electronic-structure theory density-functional theory 
(DFT) 0| is essentially limited to ground-state properties 
and its time-dependent extension [3j still displays concep- 
tual and practical difficulties. Many-body perturbation 
theory (MBPT) Q, in turn, provides a general, though 
awkward, framework for simulating electronic excitation 
processes in materials. The most elementary such pro- 
cess is the removal/addition of an electron from a system 
originally in its ground state. These processes are acces- 
sible to direct/inverse photo-emission spectroscopies and 
can be described theoretically in terms of quasi-particle 
(QP) spectra A numerically viable approach to QP 
energy levels (known as the GW approximation, GWA) 
was introduced in the 60 's but it took two decades 

for a realistic application of it to appear 0], and still now 
the routine application of MBPT to the simulation of ma- 
terials is plagued by severe numerical difficulties, which 
have limited so far these applications to models of a few 
handfuls of nonequivalent atoms, at most. The two main 
such difficulties are the necessity to calculate and manip- 
ulate large matrices representing the charge response of 
the system (electron polarizabilities or polarization prop- 
agators) on the one hand, and that of expressing such 
response functions in terms of slowly converging sums 
over empty one-electron states [1, 0, EU EH , on the other 
hand. In a recent paper, we have successfully addressed 
the first problem by expressing polarizability operators in 
terms of an optimally small set of basis functions 0|. In 
the present letter we address, and hopefully solve, the 
second problem by introducing a new method to cal- 
culate polarizability operators and self-energy operators, 
based on a Lanczos-chain technique, inspired by recent 
progresses in time-dependent density-functional pertur- 
bation theory [H, E3|- Our approach is first validated 
by calculating the vertical ionization energies of the ben- 
zene molecule, and its power demonstrated by address- 
ing the spectrum of a large molecule such as free-base 
tetraphenylporphyrin. 



QP energies (QPE) are eigenvalues of a Schrodinger- 
like QP equation (QPEq) for the so-called QP amplitudes 
(QPA), which is similar to the DFT Kohn-Sham equation 
with the exchange-correlation potential, V xc (r), replaced 
by the non-local, energy-dependent, and non-Hermitian 
self-energy operator, E(r, t',lj). In the GWA 0, 0| S is 
the convolution of the one-electron propagator, G, and 
of the dynamically screened interaction, W: 

/oo 
dw'G(v,v';uj')W(T,T';uj-uj'), 
-oo 

(1) 

where W = v + v ■ U ■ v, n(r, r';w) = (1 — P ■ v) 1 • P 
is the reducible polarizability, P the irreducible one, 
v(r,r') = jTzpi is the bare Coulomb interaction, and 
a dot indicates the product of two kernels, such as in 
v ■ n(r,r',ct;) = J dr"v(r, r")II(r", r'; uj). We assume 
time-reversal invariance to hold — so that wave-functions 
are real — and we work on the imaginary-frequency axis 
[T3 |: real-frequency results are then recovered upon ana- 
lytic continuation. One further approximation is the so 
called G°W° one, where the one-electron propagator is 
obtained from the QPEq using a model real and energy- 
independent self-energy, such as e.g. E° = V xc (r)5(r— r'), 
and the irreducible polarizability is calculated in the 
random-phase approximation (RPA): 



P°{r,r>-M = 4ReV M^M^M^MV) 
^ ilu - (e c - e v ) 



where ip and e are zero-th order QPAs and QPEs, and v 
and c suffixes indicate occupied and empty states, respec- 
tively. To first order in S' = T,a°w° — S°, QPEs are given 
by the equation: E„ w e„ + (E G o W o(E n )) n - (V X c)n, 
where (A) n = (ip n \A\ip n ), and quantum- mechanical op- 
erators are indicated by a caret. 

In Ref. [1], we reported on a strategy to build an op- 
timal representation for the polarizability operators, in 
terms of a reduced, yet controllable accurate, orthonor- 
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mal basis set {$ M (r)}: 

P°(T,r',iu) = ^P; v M$ /1 (r)$,(r'). (3) 



//// 



Using the RPA, Eq. ©, P° u (iuj) reads : 

P;„(tw) = -4Re£] L 

M e c — e v + ilj 

^^r C ^r / $ M (r)^(r)V' c (r)^(r , )'/'c(r , )*,(r , ). (4) 



Suppose that such a representation for P° can be found 
without any explicit reference to empty states (later we 
will show how this can be achieved) . In order to eliminate 
the sum over empty states in Eq. (J4}, we introduce the 
projector operator over the empty-state (electron) mani- 
fold, Q e = 1 — Qh, Qh being the projector onto occupied 
(hole) states. In terms of Q e , Eq. {4]) reads: 

p°M = 



ing Eq. ©, Eq. © reads: 

P; u (ilj) « -4Re (*a|(#° - C + i^)" 1 |t (3 )x 

Ta,v[iTfi,vv- (y) 

Having thus reduced the number of matrix elements of 
the resolvent of the Hamiltonian in Eq. ©, these ma- 
trix elements can be efficiently calculated by a Lanczos- 
chain algorithm , as explained in Ref . 13 1 . In a nut- 
shell, a frequency-independent chain of vectors is gen- 
erated for each t function. In the basis of each one of 
these chains, the Hamiltonian is tridiagonal, and thus 
easily inverted. A similar approach can be applied to the 
calculation of the expectation values of the self-energy. 
By expressing the reducible polarizability in the same 
representation used for the irreducible one, Eq. ([3|, 
n(r,r',tj) = J2u.v n Miy (iLj)$ / _ 1 (r)$ ly (r / ), the diagonal ma- 
trix elements of the correlation contribution E c to the 
self-energy, Eq. |T|), can be cast into the form: 



4Re^(^$ M |Q e (H°-e„ + ^)- 1 Q e |V^$,), (5) (t c (iw)) n = I ' du'Tl^iu' 



where H° is the QP Hamiltonian corresponding to S° 
and \ip v <& v ) is the vector whose coordinate representa- 
tion is (r|^„<I>„) = i)j v (y)'^ u (y). A direct approach to Eq. 
(0 would require the inversion of (H° — e v +iu>) for every 
value of the (imaginary) frequency and the application of 
the resulting inverse to N v x Np vectors, where N v and 
Np are the number of valence states and polarizability 
basis functions, respectively, thus making the resulting 
algorithm unwieldy. In order to substantially reduce the 
computational load, we proceed in two steps: we first re- 
duce the number of functions to which the inverse shifted 
Hamiltonian in Eq. ([5} has to be applied; in the second 
step, the inversion of the shifted Hamiltonian is avoided 
by a Lanczos-chain algorithm that need not to be re- 
peated for different values of the (imaginary) frequency 
shift. In the first step an approximate orthonormal ba- 
sis, {t Q (r)}, is built for the linear space spanned by the 
N v x N P vectors, {Q e \ipv$v)}'- 



(^„(«* M )|(fl- -i(a;-w'))~ 1 |^n(«^)), (8) 

where \ip n (vQ^)} is the vector whose coordinate repre- 
sentation reads: (r\ip n (v$^)) = V'n(r) / v(r, r')$ A1 (r')dr'. 
The matrix elements on the r.h.s. of Eq. ([8j) are calcu- 
lated with a similar procedure as for Eq. |(5]), where the 
set of vectors {\ip n (v$ n))} is first expanded into a suit- 
able optimal basis set, {s a (r)}: 



N s 



(r)S 



(9) 



(r|Q e |V>„$ M ) 



E 



(6) 



where T ajVli = (t a \Qe\ipv$n) and Nt is the number of 
t functions, which can be kept in general significantly 
smaller than N v x Np. Details on this procedure, based 
on a block version of the Gram-Schmidt algorithm, will 
be given elsewhere. Here suffice it to say that for each v 
index a block of Np vectors is first orthogonalized to the 
previously processed block, and then reduced by elimi- 
nating the eigenvectors of the overlap matrix correspond- 
ing to eigenvalues smaller than a given threshold [8j]. Us- 



and then by generating a Lanczos chain for each s; the 
convolution is finally calculated either by direct integra- 
tion or by fast Fourier transform. 

In Ref. [ij a reduced basis set for the polarizabil- 
ity operators was constructed by expressing the product 
functions, ip c (r)ip v (r), in Eq. (JH) in terms of localized 
Wannier-like orbitals. Although the number of empty 
states needed to achieve a good accuracy can be kept 
fairly small, still quite a few of them have to be calcu- 
lated. On the other hand, it was noted that this ba- 
sis can be kept independent on frequency (or on time). 
An optimal representation for the polarizability can be 
thus calculated by diagonalizing the irreducible polariz- 
ability operator at t = and keeping only those eigenvec- 
tors that correspond to eigenvalues larger than a certain 
threshold. One has: 

Po(r, r';t = 0) = J2 M^M^M^M?') 

cv 

= QhMQeM, (10) 
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Figure 1: Calculated vertical ionization potential of the ben- 
zene molecule (discs, left scale) and dimension of the polariz- 
ability basis (triangles, right scale) versus the q* cutoff. The 
polarizability basis sets have been constructed with energy 
cutoffs: E* = 5Ry (blue) and E* = 10 Ry (red). The lines 
are a guide to the eye. 

where Qh and Q e are real-space representations of the 
projectors onto the hole and electron manifolds, respec- 
tively. The eigenpairs of P (t = 0) can be easily calcu- 
lated by iterative diagonalization, noting that Q e (r, r') = 
<5(r — r') — Qh(r,r'). Such a procedure would lead how- 
ever to a number of eigenpairs much larger than strictly 
needed to achieve a good accuracy in the QP spectra. 
In order to keep the size of the polarizability basis man- 
ageable, we replace Q e in Eq. iflOj) with the projector 
over the manifold spanned by plane waves (PWs) up 
to a kinetic energy of E* , orthogonalized to the hole 
manifold (orthogonalized plane waves, OPWs), Q*. Let 
|G) = Q e \G) be one such OPW. In terms of the OPWs, 
the modified projector reads: 



Q%= ]T |G}(G'|xS G ^„ (11) 

|G| 2 ,|G'| 2 <_E* 

where 5gc = (G|G'). Using this approximation, a basis 
for the polarizability can be obtained from the eigenfunc- 
tions of the eigenvalue equation: 

V 

whose leading eigenpairs — corresponding to eigenvalues 
larger than a given threshold, q* — can be easily found by 
conjugate gradients or other iterative methods [IB]- We 
stress that this is a controlled approximation, which can 
be systematically improved; furthermore, it only affects 
the determination of an optimal representation for the 
polarizability operators, not their actual calculation once 
this basis has been determined. 

Our scheme has been implemented for norm- 
conserving pseudopotentials (PPs) in the gww . x module 



of the Quantum ESPRESSO distribution of electronic- 
structure codes [lf|, soon to be released under the 
GPL license, and benchmarked on the isolated benzene 
molecule [17j ■ In Fig. Q] we display the vertical ioniza- 
tion potential (VIP) calculated for the isolated benzene 
(GqHq) molecule with different values of the energy and 
eigenvalue cutoffs, E* and q*, defining the polarizability 
basis (see Eqs. [TT] and [12]) . The two series of calculations 
for E* = 5 Ry and E* = 10 Ry converge to the same VIP 
within few tens of meV . For both values of E* , a cutoff 
q* ~ 10 a.u. yields convergence within ~ 0.1 eV, which is 
our estimated residual accuracy due to the incertitudes of 
the analytical continuation procedure. The convergence 
of the VIP with respect to the size of the polarizability 
basis is slightly slower here than previously observed in 
Ref. 0], where sums over empty states where performed 
explicitly including a limited number of them. In Fig. 
[2] we display the differences between the calculated and 
experimental VIPs in benzene (E* = 10 Ry, q* = 0.035 
a.u., corresponding to ~ 2900 basis functions). On the 
same figure we also report results obtained: i) with a 
reduced polarizability basis (E* = 10 Ry, q* = 14.5. a.u. 
corresponding to ~ 500 basis functions); ii) with the ex- 
trapolation of the sum over virtual orbitals described in 
Ref. [i| ; Hi) with a polarizability basis of ~ 1500 PWs 
(corresponding to a kinetic-energy cutoff of 5 Ry) . Not 
unexpectedly, GW results are in good agreement with 
experiment. What matters here is that the present ap- 
proach is not only considerably faster, but even more ac- 
curate, than previous GW calculations that required ex- 
trapolation of slowly converging sums over empty states 
0. Moreover, for a same size of the polarizability ba- 
sis, the optimal polarizability basis method is also more 
accurate than the use of simple PW basis sets. 

We now demonstrate the potential of our method by 
considering the free-base tetraphenylporphyrin molecule 
(TPPH 2 ) (C44H30N4) [H]. We used a polarizability basis 
of 5000 elements, which has been obtained from Eq. lfT2|) 
using E* = 10 Ry and q* = 21.1 a.u. These parameters 
ensure an absolute convergence of the calculated QPEs 
of ~ 0.1 eV. We note that the analytic continuation pro- 
cedure involves incertitudes of similar size or even larger 
for the lowest lying states. The calculated ionization po- 
tential is 6.7 eV, in fair agreement with an experimental 
value of 6.4 eV [12]. The quality of our results is further 
illustrated in Fig. (3] where we compare the calculated 
valence electronic density of states with photoemission 
data from Ref. Nice agreement is achieved for the 

position of the peaks, while the agreement for the inten- 
sities is not as good, possibly due to our neglect of any 
matrix-element effect. 

In conclusion, we believe that the method presented 
here may open the way to ab-initio MBPT simulations 
of large and realistic models of molecular and nano- 
structured systems. While we feel that the proposed 
Lanczos technique may be considered as a sort of defini- 
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Figure 2: Differences between the calculated VIPs of ben- 
zene and experimental results. Red and green: present 
method using E* = 10 Ry, q* — 0.035 a.u. and E* = 10 Ry, 
q* = 14.5 a.u. respectively. Magenta: method of Ref. [§], upon 
extrapolation of the sum over virtual states. Blue: present 
method, but using a polarizability basis formed by PWs (see 
text). Experimental data are from Ref. [IB]. The lines are a 
guide to the eye. 
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Figure 3: Valence electronic density of states for the TPPH2 
molecule, red: theory (a Gaussian broadening of 0.25 eV has 
been used); green: photoemission data from Ref. [jg| |. 

tive answer to the sum-over-virtual-states problem, we 
think that there is still room for improving the construc- 
tion of an optimal polarizability basis. The extension 
of these ideas to real-frequency implementations of GW 
and other MBPT techniques, such as the Bethe-Salpeter 
equation, is possible, and indeed presently under way. 

P.U. thanks Xiaofeng Qian and Nicola Marzari for use- 
ful discussions and the latter for hospitality at MIT. 
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